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Abstract. 

Understanding non-equilibrium quantum dynamics of many-body systems is one 
of the most challenging problems in modern theoretical physics. While numerous 
approximate and exact solutions exist for systems in equilibrium, examples of non- 
equilibrium dynamics of many-body systems that allow reliable theoretical analysis, 
are few and far between. In this paper we discuss a broad class of time-dependent 
interacting systems subject to external linear and parabolic potentials, for which 
the many-body Schrodingcr equation can be solved using a scaling transformation. 
We demonstrate that scaling solutions exist for both local and nonlocal interactions 
and derive appropriate self-consistency equations. We apply this approach to several 
specific experimentally relevant examples of interacting bosons in one and two 
dimensions. As an intriguing result we find that weakly and strongly interacting Bose- 
gases expanding from a parabolic trap can exhibit very similar dynamics. 
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1. Introduction 

Understanding time evolution of complex quantum systems, often in the presence 
of strong correlations between constituent particles, is crucial for solving many 
fundamental problems in physics, from expansion of the early universe, to heavy ion 
collisions, to pump and probe experiments in solids. New questions of dynamical 
evolution arise in recently realized artificial quantum many-body systems, such 
as ultracold atoms in optical potentials or photons in media with strong optical 
nonlinearities. These systems are only weakly coupled to external heat baths and have 
a limited life-time, thus many experiments require interpretation in terms of coherent 
quantum dynamics rather than properties of equilibrium states. On the positive side, 
these systems allow remarkable control of parameters and open exciting opportunities 
for doing controlled experiments exploring non-equilibrium many-body dynamics. 

In the realm of many-body physics low- dimensional systems have a special place. 
They have dramatically enhanced quantum and thermal fluctuations and exhibit most 
surprising manifestations of strong correlations. Rigorous theorems provide strong 
constraints on long- range order and often such systems cannot be analyzed using 
mean-field approaches even at zero temperature. Nevertheless, equilibrium properties 
are well understood using methods specific to low dimensions, such as Coulomb-gas 
representation of vortices in two dimensions or effective low energy descriptions of one- 
dimensional systems including Luttinger liquid and sine-Gordon models (see e.g. ref. 
PQ). However, such analysis cannot be straightforwardly extended to non-equilibrium 
dynamics. Most equilibrium theories focus on the low-energy part of the spectrum 
while non-equilibrium dynamics can couple degrees of freedom at very different energy 
scales ESIEISIIHIIHIITIIHIIHI [101 [TH [121 [131 H31 US]. It would be highly valuable to 
have examples of many-body dynamics of low- dimensional strongly correlated systems 
amenable to an unbiased analytical treatment. These examples could be used not 
only for analyzing experimental systems, but also for testing theoretical calculations 
utilizing effective models or approximations and for checking validity of new numerical 
approaches. In this paper we propose such a class of non-equilibrium quantum problems 
with time-dependent Hamiltonians which allow for a scaling ansatz of many-body wave 
functions. 

Scaling solutions in quantum dynamics were first discussed in the context of a 
single harmonic oscillator with a time-dependent frequency [161 ETJ EES EEH1 EQl |2"T] . 
This problem can be reduced to a time-independent one by properly rescaling space 
and time. Scaling transformation of variables is possible due to the existence of a 
dynamical symmetry generated by dynamical invariants of the system [181 EH] • There 
are also extensions of this approach to single particle problems with potentials of the 
Coulomb and inverse square type [221 E31 [17]. In the context of many-body problems, 
scaling has first been used within mean-field approaches to bosonic systems, for the 
classical Gross-Pitaevskii equation pH |25l EHl [29l [30] . Beyond these effective one- 
body problems, scaling solutions exist for hard-core bosons in one dimension [31] and 
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in the unitary limit of fermionic gases with infinite scattering length these are 
problems for which the interaction enters a constraint on the wave function of an 
otherwise non-interacting system analysis. Away from these specific limits, Pitaevskii 
and Rosch [27J introduced a scaling ansatz for a two-dimensional many-body system of 
particles interacting with contact or inverse square interaction and related the existence 
of such solution to a hidden SO(2, 1) symmetry. In this paper we further extend 
full many-body scaling solutions to more general types of interaction and arbitrary 
dimensionality. This generalization can be achieved by allowing three parameters of the 
system - the mass, the interaction constant and the external potential - to be time- 
dependent. Scaling solution is possible when the interdependence of these parameters is 
given by an Ermakov type equation, similar to the one discussed in earlier approaches 
[3T1 l26l [27] . and an additional self-consistency equation which depends on dimensionality 
of the system and the nature of interactions. 

Dynamical control over the system parameters is possible in recently developed 
artificial quantum systems, such as trapped ultracold atomic gases, where the effective 
interaction can be tuned using either Feshbach resonances or by changing the transverse 
confining potential, whereas the effective mass can be changed by application of the weak 
optical lattice [32]. Also with photons in nonlinear optical devices, where the time- 
dependent dispersion and Kerr nonlinearity can be achieved using electromagnetically 
induced transparency [331 EH EHJ EHJ [37]. In this paper we propose applications of the 
scaling ansatz which are experimentally relevant in the context of both of these systems. 

We emphasize that apart from the tunability of the parameters no specific 
restrictions on the system properties are imposed. Particles can obey fermionic, 
bosonic or mixed statistics, interact by pairwise interaction, and be subject to parabolic 
confining potential, to a linear potential, and to a complex chemical potential. The 
basic idea of the scaling solution presented hereafter is to map the non-equilibrium 
equations of motion to an equilibrium many-body Schrodinger equation. The mapping is 
based on scaling functions which relate correlation functions of time-dependent systems 
to correlation functions of systems in equilibrium. Hence several results known for 
equilibrium many-body systems can be directly translated to non-equilibrium situations. 
The reverse conclusion is also true: from a measurement of the system out of equilibrium, 
e.g. a quantum gas after expansion, we can deduce its initial (equilibrium) properties 
[38]. 

The paper is organized as follows. In section [2] we introduce a general formalism 
of scaling transformation for a many-body Schrodinger equation. In section [3] as an 
example of application of our approach we compute momentum distributions for one- 
and two-dimensional bosonic gases with contact interactions released from a parabolic 
trap. Further details are given in the Appendices, where we also discuss relation of 
our work to classical integrability of time-dependent bosonic systems with contact 
interactions. 
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2. Scaling transformation — general approach 

Our starting point is the many-body Schrodinger equation for N interacting particles 
in D dimensions, 

— = H(t)V(x 1 ,...,x N ;t), (1) 

i N N 

1=1 l^j 

where Ajj?) is a D- dimensional Laplacian acting on the coordinate x, = 
{xp ,Xi*\ ■ ■ ■ ) of the particle i (h = 1 here). The external parameters (chemical 
potential f/,(t), linear potential g(t) and trapping frequency oo(t)) and the many- 
body interaction potential V^(x; t) depend explicitly on time. The chemical potential 
u(t) = 9? [//(£)] + i^s[fi(t)} can accommodate effects of dissipation via its imaginary part 
While the dependencies on the linear and chemical potentials can be removed by the 
Gallilei transformations and phase shifts respectively, we note that solving the quantum 
problem with time dependence of the remaining parameters represents a non-trivial 
task. For instances, unlike in the non-interacting case, the time dependence of the mass 
can not be removed by the simple redefinition of time variable. 

We address the following question: under which conditions Eq. (JTJ (the ^-system) 
can be transformed into the Schrodinger equation for a time-independent (<£-) system: 

,a$(y 1 ,...,y JV ;r) 



dr 



H ^( yi ,...,y N ;r), (2) 



i=l i 

We emphasize that so far in (j2J) Wo an d tuq are unspecified parameters; in particular the 
$-system can have vanishing confining potential even when the \l/-system is confined. We 
assume that the time dependence of the pairwise interaction potential enters through a 
single time-dependent coupling V(x;t) = V(x.)v(t) and Vo(x) = V(yi)v . We further 
assume that the interactions have a scaling property and are characterized by the 
exponent a, which we take to be the same for both \&- and ^-systems, 

V{Xx) = \ a V(x). (3) 

Most generic interaction potentials (or pseudo-potentials) satisfy a scaling law ([3]): s- 
wave interactions V s {x) oc 5{x) (a = —D), any algebraic law, V(x) ex |x| a , including 
Coulomb (a = —1), inverse square law (a = — 2) or dipole-dipole interactions (a = —3). 
Other examples are ultracold fermions interacting via p-wave channel which gives rise to 

| It is known that the time evolution under non-Hermitian Hamiltonian in a spirit of stochastic wave 
function description is equivalent to the description of the open system by the Lindblad master equation, 
see e.g. Ref. [37] 
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the 5' pseudo-potential (a = D — 1). Also logarithmic potentials can be treated; scaling 
of the logarithmic law produces a time-dependent shift to nit). 

To express the solution of the time- dependent Schrodinger equation ([T]) in terms of 
the solution $(yi, . . . , y^; r) of the static equation fl2]) we introduce the scaling ansatz 

tt( Xl , . . . , Xjv ;t) = e WEf =1 *?+G( t )£f =1 x i+ M( t) ;v] 

x ^^y*(yi,---,yjv^), (4) 



with = (xi/L(t)) + S(t) and r = r(t). Direct calculation shows (see Appendix A ), 
that this ansatz is valid if the scaling functions R(t), L(t), F(t),r(t),G(t),S(t), M(t) 
satisfy a set of coupled differential equations, 

R(t)= -J--DF(t)R(t)-^(t)}R(t), (5) 

L(t) = -^-F(t)L(t) (6) 
m(t) 

F{t) ~ ~ MT) F {t) 2 + 2L*(t)m(t)' (?) 

M(t) " " 2^5) " W)] + 2m(t)L 2 (t)' (9) 

6(*) = -2?^ + (ii) 

m(£) m{t)L 6 {t) 

L -^\t) = (i 2 ) 

m u 



It is not obvious a priori that equations (151 112]) can be satisfied simultaneously for any 
reasonable time-dependencies of system parameters m(t), v(t), u{t). Our next goal 
is to show that there is a number of non-trivial cases for which equations fl5l fT2|) are 
consistent with each other. First of all we note that equations (jSJ) and (JH]) imply that 
R(t) = [L(t)p/ 2 exp(- J * $f\pi(t)]dt). In the absence of dissipation (3 [//(*)] = 0) this 
condition is equivalent to the conservation of the norm of the wave function under the 
scaling transformation. Eq. ([6]) allows to express F(t) via Lit), Fit) = ?p-L/L, which 
can the be substituted into the Eq. ([2l). This leads to the differential equation for L(t), 

L[t) + h(t)L(t) + u\t)Lit) = ^_|_ , (13) 

where h(t) = m m(£)/m(t). The term with the first derivative can be removed by the 
change of variables Lit) = exp[B(t)]y(t) with Bit) = —h/2. For y{t) we obtain 

yit) + ^(t)yit) = ^f- (14) 
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where Q 2 (t) = |/i 2 — |/i + u 2 (t). Eq. (fT4"l) is the celebrated Ermakov equation [39] first 
discovered in 1880 140] . This equation has been used primarily for tracking invariants of 



the time-dependent harmonic oscillator. In Appendix B we show how one can use the 



non-linear superposition principle to reduce Eq. (fl4|) to the linear equation. Once L(t) 
is known, the remaining set of equations for S(t), M(t), G(t) can be solved directly. 

In summary, to find time-dependent parameters which admit a scaling solution 
one can apply the following recipe: after specifying two time-dependent functions 
u(t) and m(t) one obtains a solution of the Ermakov equation ( fT4j) from which one 
determines time-dependent interaction strength v (t) consistent with Eq. (fT2|) . Solutions 
for the functions M(t), G(t), S(t) can then be obtained straightforwardly provided that 
functions g(t) and fi(t) are explicitly specified. Note that the complexity of our method 
(e.g. solving the Ermakov Eq.) does not depend on the number of particles N. 

The initial conditions for systems (OQ) and (j2J) are related to each other through 
Eq. applied at time t = 0: 

*( Xl , . . . , XA r; 0) = e *[^(o)Et 1 x?+G(o)E l i I 1 x l+ M ( o)iv] 

R N (0) vL(0) ' • • • ' L(0) ' K } )' 

Generally at t = the Hamiltonians controlling the dynamics of and ^-systems do 
not coincide. For example, they can have different confining potentials, or one system 
can be in a trap while the other one is in free space (oo = 0). In this paper we 
focus on a finite initial trapping potential, oj(0) = ojq > 0, for which we introduce the 
additional assumption that at t — the two systems coincide. This means that we 
have m(t = 0) = m , v{t = 0) = v F(t = 0) = G{t = 0) = M(t = 0) = 0. At 
t > the parameters of the ^/-system begin to change in time while the parameters 
of the ^-system remain constant. Since the two systems coincide for t < 0, the initial 
state of the \l/-systems at t — should correspond to the equilibrium state of the $- 
system. Existence of the scaling solution in one dimension in the hard-core limit v — > oo 
has been established previously [3T]. Within our approach this can be understood as 
follows: the first equation of (|T2|) is trivially satisfied, whereas other equations do not 
depend on the interaction strength and remain valid. Another special case is the two- 
dimensional system with contact interactions studied previously by Pitaevskii and Rosch 
(D = —a = 2), for which Eq. (JT2J) is satisfied by constant mass and interaction. 



3. Dynamics of Bose-gas with contact interaction released from the trap 

In this section, as an example, we apply the scaling approach to an ultracold Bose gas 
with contact interaction which is prepared in a confined, weakly interacting initial state. 
The nontrivial dynamics comes from a sudden switching off of the confining potential 
from u)(t) = uo at t = to u(t) = at t > 0. Solution of the scaling equation (113"]) 
for constant mass m(t) = m , is then given by L{t) = \J (1 + ^ot 2 ), and consequently 

F(t) = ^f^/L 2 (t). In |Appendix B~3 we examine additional scenarios corresponding 
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to varying mass which exhibit similar behavior of the scaling functions. Here, we also 
assume that /i(t),g(t) are time-independent constants. 

To characterize the non-equilibrium dynamics it is convenient to deal with 
correlation functions which can be easily derived within the scaling approach 
( Appendix D ). The dynamics of the momentum distribution, for example, can be related 
to the single-particle density matrix gi of the initial state, 

/oo roo 
dx / dx / 0i(x,x';O) 
-oo J —oo 

x e -i[F(f)L2(t)(x=-x' 2 ) + l(t)p( X -x')] _ ( 15 ) 

From the asymptotic behavior of the scaling functions L(t) — > WQ t>iWo^ and F(t)L(t) = 
m(t)L(t)/2 — > wot>1 mou;o/2 we can extract the long-time limit of the momentum 
distribution using the stationary phase approximation (SPA), 

n(p,t)-> f — ^ — ^ — , ^ — ;0). 

> «*,*>! \ m (t)L{t)J m{t)L{t) m{t)L(t) 

Hence the momentum distribution becomes fully determined by the density distribution 
[p(x, t) = #i(x, x, t)] of the initial state. 

For a quantitative description of dynamics we need to specify the initial correlation 
function, which we take from earlier analysis of effective theories for weakly interacting 
Bose-gases in harmonic traps [H I l4"2 l 14*3]. An important characteristic for a condensed 
state with a sufficiently large number of particles is the Thomas-Fermi shape of the 



density profile, p(x) = Q(Rtf — l x l)(/ i /^o) (1 — ( x / Rtf) 2 ), where Rtf = y 2/V m oA^o 
is the Thomas-Fermi radius. 

First we analyze the one- dimensional case in the low-temperature regime when the 
coherence length is of the order of the Thomas- Fermi radius (Eq. (1E.1|) of Appendix D ). 
According to the scaling equation ( [12]) , for contact interactions, ^(x, t) = v(t)5(x.) 
(a = —D), the interaction must be tuned inversely proportional to the scaling function, 
v(t) = Vq/ L{t). In Fig. [T^i results of numerical evaluation of the momentum distributions 
( fT5|) for specific initial values are shown together with results from SPA. The behavior 
of the p = component is characterized by a steep decay on a time scale Uq 1 followed 
by slowly dephasing oscillations, which are due to the finite extension of the density 
profile and the quadratic phase factor in ( fl~5l) . The corresponding period of oscillations 
P is determined by the Thomas- Fermi radius, P ~ fpr 2 "- Oscillations as a function of 
|p| at constant t can be attributed to the finite Thomas- Fermi radius as well. Here the 
quadratic phase factor leads to the oscillation period growing with |p|. In agreement 
with the SPA prediction, the momentum distribution relaxes to a semi-circle law. This is 
remarkable, since such a behavior has been previously associated with one-dimensional 
Bose-systems in the strongly interacting limit (v — > oo) p]T] only. In our case the 
interaction strength is initially small and then even decreases in time. We note that this 
can not be understood as effect of dilution due to expansion of the system because the 
effective one-dimensional interaction parameter [44], 7 oc v(t)/p(t) oc v(t)L(t), remains 
constant. 
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Figure 1. Temporal evolution of momentum distribution functions following turning 
off the trap at t = 0. The insets show the time evolution of the p = component. The 
initial correlation functions are derived from effective theories (Refs. [4ll [42j [43] , see 
also |Appendix D[ ) . Dynamical evolution is obtained from numerical integration of Eq. 
(|T5|) . The stationary phase approximation (SPA) represents the asymptotic i ^ oo 
result. Numerical errors are of the order of the line thickness. In the one-di mensiona l 
case (a) the system parameters are N = 140, fc^T = Q.IHloq, vq = 0.2 fi 3 wo/mo, 

Rtf = 3.46y^/ (mowo), v(t) = vq\/(1 + uj^t 2 ). In the two-dimensional case (b) the 
interaction strength is constant, v (t) = vq and N = 16, ksT = O.ltiajQ, vq = 0.27i/mo, 
Rtf = lAly/h/(m ujo). 



In two dimensions a = —2 and Eq. ( Tl~2"l) leads to interactions which are constant in 
time. When the initial state is weakly interacting (Appendix D), we choose an effective 
theory which incorporates effects of quantum and thermal fluctuations [12]. Results of 
numerical evaluation of Eq. (TH)]) are shown in Fig. ([Do)- The momentum distribution 
evolves very much like in the one- dimensional case and is essentially determined by the 
initial density distribution and the associated Thomas-Fermi radius. Here the number 
of particles (N = 16) is set to be smaller than in the one-dimensional system. Therefore 
the asymptotic stationary phase solution is approached slowly and oscillations dominate 
in the analyzed time window uot < 20. We checked that both in one and two dimensions 
the results are robust against variation of temperature and interactions as long as phase 
coherence is not destroyed. 

The analysis of these examples leads to remarkable consequences. We note that 
the stationary phase regime is reached rather quickly with momentum distribution 
determined by the initial density distribution. Therefore specially designed initial 
density distributions (equilibrium or not) can be used to create specific momentum 
distributions, such as step-like fermionic ones, on demand. It is remarkable that such 
behaviour, which has been obtained previously in the strongly interacting limit, persists 
down to arbitrarily weak strength of interaction. This is opposite to what is realized 
in time-of-flight experiments of ultracold atoms released from a lattice jl5], where the 
expansion at sufficiently large times can be regarded as free and momentum distributions 
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get mapped to density profiles. By contrast in our case we find that the real space density 
profile in the trap determines momentum distribution after expansion (see Eq. (fT5|) ). 
While we do not discuss the appropriate time evolution of u(t), m(t), and v(t) here, we 
point out that the time-of-flight 'far-field' limit [15] may also be captured formally by 
our scaling approach when the asymptotics of L(t) are linear and the contribution of 
the quadratic phase factor in Eq. ffT5l) . m(t)L(t), vanishes in the long-time limit. 

4. Conclusions and outlook 

We used scaling ansatz to show that certain quantum non-equilibrium problems with 
time-dependent parameters can be related to equilibrium problems with constant 
parameters provided that the time-dependent parameters satisfy a system of self- 
consistency equations. This approach is valid for rather general types of interactions 
and is not linked to the integrability of the model. However, an integrable structure, 
when it exists, is consistent with the scaling transformation. Solvability by the scaling 
ansatz is a consequence of the non-relativistic dynamical symmetry which received 
considerable attention recently in relation to the non-relativistic version of AdS/CFT 
correspondence [IHl E7J HBJ EHJ ED]. The appearance of this symmetry in realistic many- 
body systems, which we discuss in this paper, can open intriguing connections to the 
concept of AdS/CFT correspondence. 

We used scaling approach to analyze the problem of an abrupt switching off of a 
confining potential for bosonic systems with contact interactions in d — 1 and 2. Such 
experiments can be performed using either ultracold atoms or photons in non-linear 
medium. We find that the asymptotic momentum distribution is essentially given by 
the initial density profile - a phenomenon which previously has been discussed only in 
the (Tonks-Girardeau) limit of the infinitely strong repulsive one-dimensional Bose gas 
|31j . Possible future applications of the scaling ansatz include interaction quenches or 
transport phenomena (by considering finite linear potentials). Extensions of our method 
to systems with dissipation are also possible. 

In our analysis we considered the situation when the scaling ansatz is obeyed 
exactly. We expect however that our results remain qualitatively valid even for systems 
with small deviations from the exactly scalable Hamiltonians. For example, weak lattice 
potentials should not have dramatic effects as long as the effective mass approximation is 
applicable. Therefore one could achieve a full description of time-of-flight experiments if 
the lattice potential and interactions are tuned accordingly. Moreover it is conceivable 
that on a phenomenological level the ansatz can be used even when the time- and 
space-dependencies of system parameters do not fully satisfy the consistency equations. 
The scaling solution could then be seen as a universality class of non-equilibrium 
systems, very much like a renormalization group fixed point at equilibrium. It would be 
interesting to address this conjecture in experiments. 
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Appendix A. Derivation of the scaling equations 

We consider the ansatz (jlj) 

^ N N 

*(x l5 . . . , x^; t) = — - exp(z[F(t) ]T xf + G(t) £ x, + M(t)]) 
K W i=i i=i 

for the transformation between the many-body Schrodinger equation with time- 
dependent parameters (Eq. (jTJ and the equation (J2J) with time-independent coefficients. 
Calculating directly 

<p=( -^ + f^ x ' +! f^ x,+ 4 y*""'^- t » < A - 2 > 

+ i e *-> f; lx ,(- A) + s (i) ] + y^^^f , 

where for the sake of brevity we introduced </>(xj, t) = F(t) Y^iLi x i +G(£) X^i X i + M(t) 
and where the dot denotes the derivative with respect to t, and 
0*(xi,t) 1 



chc; R 



2tf^x ! + G e*'^(y„f) (A.3) 



, 1 j^ dHyur) 
AgpVfat) = { + ^(2zFx, + zG)(2zF Xi + iG)) $(y i; r) (A.4) 



r , riG ^ ;t) +Ag)<t(yt;() M }e 



itj)(xi,t) 



\ RL dy, Vi w ' ' RL 2 

Substituting this into the initial Schrodinger equation (CQ) with time- dependent 
coefficients and adding and subtracting the term A(t) Si x ? with yet to be determined 
function A(t) we regroup the different contributions in front of $(yj,r), <9$(yj, t)/<9vj, 
and A y .. Each group has several contributions proportional to x°,Xj,x^ which are 
linearly independent and must be treated separately. This is how conditions expressed 
by Eqs.((7j) appear. The remaining equation has the form of a Schrodinger equation with 
time-dependent coefficients 

M{yi,r). l 



% ^ T 



dr 2m(t)L 2 (t) 



A w $(2/i,r) (A.5) 



+ 



A(t)L 2 (t)Y,yl + L a (t)v(t)V(y t - yj : 
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We note that to compensate the terms appearing after the change x, — > y, in the 
quadratic potential we get terms proportional to Uq in the Eqs. fl6l fT2|) . Now, requiring 
that the three unknown functions r, L(t),A(t) satisfy 

f = UW¥y W = <*)V®> A(t)L\t) = r^fl (A.6) 



we obtain the remaining conditions in the set of Eqs.(!5T ll2]) . Under this conditions the 
Schrodinger equation for the function $(?/, r) has no time- dependent coefficients. From 
the conditions (1A.6|) above we determine the function 



4 
a + 2 



A(t) = ^<4Mt)rnW (A.7) 
2m(t)vo +2 

Therefore we find that when pairwise potentials obey Eq. ([3]), and the systems of Eqs 
(ED-flE! is satisfied, Eq. flD is indeed mapped to Eq. Q2J). 

Appendix B. Analysis of the scaling equations and their solutions — the 
Ermakov equation and dynamical symmetry 

Appendix B.l. General properties of the Ermakov and related equations 

In this Appendix we briefly overview some general properties of the Ermakov (sometimes 
spelled as Yermakov) equation which plays such a fundamental role in our formalism. 
We also point out the relation of this equation with the Riccati equation and with 
the linear differential equation with variable coefficients. The Riccati equation directly 
appears in our approach in some limiting cases. 

The Ermakov jlQ] equation is defined as follows 

y{t) + f(t)y(t) = (B.l) 

Here a is some t- independent constant. If there is a nontrivial solution of the second 
order differential equation 

x{t) + f{t)x{t) = (B.2) 

then the transformation 

puts the Ermakov equation into the form 

z# = az~ 3 . (B.4) 

where the subscript denotes the derivative. The solution for the initial equation then 
follows immediately 



dt. 



CV = ax 2 + x\C 2 + d / (B.5) 
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where are arbitrary constants. If we take two solutions of the linear (Hill) equation 
to satisfy initial data 2i(0) = x±, ii(0) = ii while x 2 (0) = 0, x 2 7^ then a general 
solution of the Ermakov equation is given by a nonlinear superposition principle, 



y(t) = Jx\{t) + — 9 xl{t) (B.6) 



w 2 



where w = X\Xi — x<iX\ is a constant Wronskian. 

Now, provided the linear equation for x(t) is satisfied, the function u(t) defined as 

x(t) = exp(- f u{t)dt) (B.7) 
Jo 

satisfies the Riccati equation, 

u-u 2 = f(t) (B.8) 

This demonstrates that all three equations are closely related: Ermakov, linear 
second order differential equation with variable coefficients and the Riccati equation. 
Other remarkable equations are also connected to the Ermakov equation. For example 
(taking a = 1 for simplicity in ( 1B.4I) ) and defining ^(t) = z(t)~ 2 we obtain ££— (3/2)(£) 2 + 
2£ 4 = 0. Now, defining w(t) via £(t) = aw/w with a 2 = —1/4 we obtain a Kummer- 
Schwarz equation ww — (3/2) (iv) 2 = 0. 

In some limiting situations (e.g. ojq = 0, see the next appendices) the Riccati 
equation appears naturally in our approach, so we sketch some of its properties here. 
The general Riccati equation with time- dependent coefficients 

ii(t) = f(t)u\t)+g(t)u(t) + h(t) (B.9) 

can be transformed into the second order differential equation 

f(t)y(t) - [f(t) + f(t)g(t)]y(t) + f(t)h(t)y(t) = (B.10) 

by the following substitution y(t) = exp(— / f(t)u(t)dt). In many cases a particular 
solution of flB.101) is easier to find than the one for the flB.9j) . 



The Riccati equation has a remarkable property: if there is a known particular 
solution u (t) of (IB. 91) . then the general solution of (IB. 91) is given by 



u(t) = u (t) + $(t) 



c-j f(t)$(t)dt] 1 



$(t) = exp 



(2f(t)u (t)+g(t))dt 



(B.11) 
(B.12) 



where C is an arbitrary constant. The particular solution uq(x) corresponds to C = 00. 

The property (IB. Ill) allows the construction of many solutions of (IB. 91) for given 
functions / (t) , g (t) , h(t) . If, for example, f(t) = 1, g(t) is arbitrary and h(t) = 
— (a 2 + ag(t)) a particular solution is %(£) = a , an d a general solution is then 

u{t) = a + $(t)[C - J ^{t)}-\ $(t) = exp(2at + J g{t)dx) (B.13) 

for arbitrary C. For example for f(x) = 1, g(x) = 0, h(x) = bx n we obtain 

u{t) = w(t) = ^[d Jx (iVtt*) + C 2 Y, {\Vbt% (B.14) 

w(t) 2k k 2k k 
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k=-(n + 2), for n^2 (B.15) 

u{t] ) = T~ t2X{ ; 2ATT t2X + c ^^ for n = - 2 ' ( R16 ) 

where A is a root of A 2 + A + b = 0. 

Appendix B.2. Relation to dynamical symmetry 

The Ermakov equation has the symmetry algebra isomorphic to sl(2,R), which is 
isomorphic to the algebra so(2, 1) of rotations on the surface of one-sheet hyperboloid. 
The property f jB.llj) of the Riccati equation is related to the covariance of the Riccati 
equation with respect to the fractional- linear transformations which are generated by 
the action of sl(2,R) algebra: the general solution can be expressed as a combination 
of particular solutions. The same algebra (more explicitly, one of its form, su(l, 1)) 
appears as a dynamical symmetry of the quantum harmonic oscillator, where Ermakov 
equation appears as well. This has been first found in [16J. There a single quantum 
harmonic oscillator with time- dependent frequency has been solved using the methods 
of (adiabatic) invariants. An adiabatic invariant in this case is a function of a solution 
of the Ermakov equation. This approach has led to appearance of the Ermakov- Pinney 
type equation [40] in quantum mechanics (see e.g. [20] for a recent review). In [17] 
the same equation appears as a certain consistency condition on the time-dependent 
rescaling of coordinate and time in the wave function of the oscillator. It became 
clear that these two approaches, one based on dynamical invariants and the other 
on the scaling of dynamical variables, are equivalent. Indeed the rescaling procedure 
can be regarded as a transformation, generated by a certain symmetry group, i.e. 
s/(2,i£). The generators of this symmetry are operators corresponding to dynamical 
invariants. Therefore the successiveness of applicability of scaling transformation implies 
the presence of dynamical symmetry generated by the dynamical invariants [TBI 119] . 
For this symmetry to hold one has to have a special class of potential terms in the 
single-particle Hamiltonian [22]. Physically interesting potentials correspond to the 
contact interaction, harmonic, Coulomb and inverse square laws. That is why the 
scaling approach has been applied to a Calogero- Sutherland model [23] and classical 
Gross-Pitaevski type systems [2U [251 [23 [281 ED]- The appearance of the su(l, 1) 
dynamical symmetry in our non-relativistic systems suggests a possible connection to 
non-relativistic version of the AdS/CFT correspondence [16] , [UJ , [H] , [H] , [50] . In fact 
the Virasoro algebra of any conformal field theory contains su(l, 1) as subalgebra. 

Appendix B.3. Specific solutions for uq > 

We compare examples for decreasing trapping potential and constant, increasing and 
decreasing masses. 

(a) Constant mass - For the case of constant mass m(t) = mo we choose an exponential 
decrease of the potential co(t) = Woe - '^". The two independent solutions of the 
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Figure Bl. Scaling functions for ujq > (h reinserted by dimensional analysis). Each 
curve corresponds to one of the cases (a)-(c) analyzed in the text. 



homogeneous equation flB.21) read Xi(t) = J (2r w y u>(t)), x 2 (t) = Yq{2t^ujqJ oj(t)). 
In fig. IB 1 1 the resulting scaling functions obeying the initial conditions L(0) = 1, 
F(Q) = are plotted. For sufficiently small r u the functions are well described by 
the limit r w — > 0, for which the scaling solution reduces to 

L(t) = V / (T7^), F(t) = ^/L\t) . (B.17) 

(b) Increasing mass - We choose m(t) = m§e t l Tm and, for sake of simplicity, tu(t > 0) = 
0. The solution then reads 

Lit) = ^\ + {l-e^)rlm ul 

F(t) = (1 - e-^rl m u 2 /L 2 (t) , (B.18) 

(plotted in fig. IB1I) ; this is similar to the scaling functions of the case (a), although 
the time is rescaled and in the limit t — > oo the functions converge to the values of 
the functions of case (a) at t — T m . 

(c) Decreasing mass - For m(t) = m§e~ l l Tm the scaling functions take the form of case 
(b) when replacing r m by — r m (see fig. IB II for an illustration). 

We emphasize that the solutions do not depend on the dimensionality of the system; 
only the interaction constants, which have to fulfill the consistency equation ( jl~2l) . will 
do so. 

i 



Appendix B.4- Specific solutions for ojq = 

Based on two examples we demonstrate within our formalism, that if we relate the non- 
equilibrium system in the trap to the system without trap (the case Uq = in the main 
text) we directly obtain a Riccati equation. 
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For D = 1 Eq. [Threads Lit) = m (m(t)c(t)) _1 (we define c(t) = v(t)/vo) what 
we substitute in the equation for L{t) to obtain F(t) = —(m(t)/2)jt\og[c(t)m(t)/mo\ 
. Consistency with the equation for F(t) imposes the following relation between three 
time-dependent parameters 

rh(t) d , ... . m(t) d 2 , r . , . < > 

- -^ Jt \og[c(t)m(t)\ - -Y—\ g[c{t)m{t)] (B.19) 

?n(t) . d , r / \ / vni m(t)uj(t) ,„ 
= ~ 2 ( aY I °8[ c (*)»n(*)]) a - 2 • (B.20) 

By introducing C/(i) = 4 log(c(t)m(t)) it reduces to the Riccati equation 

17(f) = w(f) - -^(log[m(t)])C/ + f/ 2 . (B.21) 
at 

The scaling ansatz (j3J) implies the relation between initial conditions of the two systems: 
= 0) = exp(zF(0) Ei^ 2 )*(^ = 0) provided that L(t = 0) = 1. The initial condition 

for the function U (t) is not so important for us because of the special property of the 

Riccati equation, related to the Backlund symmetry, which allows to interrelate solutions 

with different initial conditions via a rational function. 

We note that the same equation describes the evolution of spin in a time-dependent 

magnetic field. A general way to solve it is to notice that under some change of variables 

it can be reduced to the second-order liner differential equation 

u - P[t)u + Q{t)u = 0, P(t) = -j t log[m(f)], Q(t) = u(t). (B.22) 

Numerous explicit solutions are possible if we specify the functions u(t),m(t). 

In the two-dimensional case we obtain from (!5])- (jl2p that R(t) = L(t) and 
time-dependent parameters are connected by the constraint c(t)m(t) = Cq. Then 
F(t) = (m(0/2)|log[L(t)]. Introducing V(t) = | logL(t)) and h(t) = |log(m(t)/2) 
we obtain 

-^^ = u(t) + h(t)V(t)+V 2 (t) (B.23) 
dt 

which is a Riccati equation for the coordinate scaling function L(t); its solution for 
given time-dependent parameters m(t),tu(t) then defines a solution for the time-rescaling 
function 

dr(t) m 



dt " m(t)L 2 (t) 
To be specific we list two examples of dynamical parameters: 



(B.24) 



Increasing mass - From the form of the Riccati equation it is somewhat appealing 
to take m(t) = moe at , and constant u(t) = Q. Then 

c(t)=<j>(t)exp[-acb/2] (B.25) 

where <f>{£) = sin(At + B)/C with A, B,C related to a and Q. In particular, for 
m{t) = e 2t , where Q = 1,A = B = C and C ->• we obtain c(t) = (1 + t)e _i . 
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(b) Constant mass - For m(t) = m the equation can be transformed into the equation 
for the harmonic oscillator with time-dependent-frequency oj(t) for which many 
known solutions exist. Using these solutions we can extract the function c(t). In 
particular, for constant u(t) = Q the solution for some domain of parameters is 

c(i) = — — . (B.26) 

v ' m cos(fit) 

In the simplest case of m(t) = 1, cu(t) = we obtain c{t) = —1/(1 + t). This 
example is a many-body analogue of the solution of the Hamiltonian with potential 
V(x) = c(t)5(x) found in ref. [22] for a single-particle Schrodinger equation. Direct 
application of this solution can be found in the ultracold Bose gas close to the 
confinement- induced resonance I5T1. 



Other examples of solutions of (1B.9|) can be found in the literature, see e.g. ref. 



1521. 



Appendix C. Classical integrability of the nonlinear Schrodinger equation 
with time-dependent parameters 

It is instructive to check whether the exact scaling transformation we have studied in 
this paper is consistent with the property of integrability of the nonlinear Schrodinger 
equation (NSE). Here we address this question for the classical NSE. 
In the zero curvature representation, the NSE 

d 2 ^ 

i-^—^cW* (C.l) 
is represented by the system of the first order differential equations 



^ = U(x,t,X)F, ^ = V(x,t,X)F, F =( f f l) (C-2) 

such that the matrices U (x, t, A) and V(x, t, A) which depend on the spectral parameter 
A satisfy the condition 
dU dV 

which is equivalent to the compatibility condition of the system, 

— = — (C 4) 

dxdt dtdx 

and which is equivalent to the initial Schrodinger equation. In case of (10. lft one can 
establish that 

U =U + \U U V = V + XVi + \ 2 V 2 (C.5) 

U = yfc{^a + + *tr_), Ui = ^ia 3 (C.6) 

V =tc\<S>\ 2 o- 3 -z^- C (—o- + -—o-_), V 1 = -U , V 2 = -U 1 (C.7) 



Scaling approach to quantum non- equilibrium dynamics of many-body systems 



17 



Conserved quantities are constructed from the matrices U, V in a known way. This 
method provides a direct way to various generalizations of NSE. In particular one can 
obtain some generalization where the interaction parameter c and the mass are explicitly 
time-dependent functions. Introducing generalization of (1C.5|) as 

f -\a{x,t) -Y(x t t)V \ 



U 



-y(x,t)V POM) 



z^(|M>| 2 , A(x,t)) B{^,^,^x,t)) 



" ^ B*(*,f -^(i^i 2 ,A(x,t) ; 1 • } 

one can look for generalizations of integrable NSE by appropriately choosing the 
functions a(x, t), /3(x, t), ^(x, i), X(x, t), ^(x, t), A, B, D. Analysis of the zero-curvature 
condition (lC.3j) in the case of inhomogeneous time- dependent functions leads to a set of 
equations between those functions and reveals a large class of solutions of the classical 
equations of motions for NSE with time-dependent coefficients. To get a consistency 
condition for a zero-curvature representation we conclude that the spectral parameter 
should be an inhomogeneous time- dependent function. 

Some restricted form of this inhomogeneous time-dependent U — V pair has 
been considered in ref. [53] where it was shown that a combination of space-time 
transformation together with a U(l) gauge transformation of the linear equations for the 
U—V pair and corresponding redefinition of the field variables brings the system into the 
form of a homogeneous time-independent NSE system, thus showing the integrability of 
a time-dependent system. We note that a similar analysis has been given in Ref. [51]. 

Although it is more difficult to show integrability on the quantum level directly, 
presumably the property of integrability is not violated in that case for specific choice 
of time-dependent parameters which correspond to our scaling equations. A related 
approach based on the inhomogeneity of spectral parameters for the quantum sine- 
Gordon model has been recently presented in ref. [55] . 



Appendix D. Scaling of correlation functions 

With the scaling ansatz (jl]) the relation between the single-particle correlation functions 
in the time-dependent and time-independent systems is derived straightforwardly, 

poo poo 

#J j (x,x',i)= / .../ dx 2 ...dx JV **(x,x 2 ,...,xj V ;t) 

J —oo J —oo 

x *(x', x 2 , . . . , xjv; t) 

The labels in the gi-function refer to the time-dependent (^) and time-independent 
($) systems. From this expression we can readily extract the density: p^(x, t) = 
<7i(x, x, £) = (1/L(i))p(*)(x/L(i); 0). The momentum distribution of a time-dependent 
system, defined as 

dx / dx'e- ip ( x - y )g{*\x,x',t), (D.2) 

-oo J — oo 
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is then given by 

/oo r oo 

-oo J —oo 

exp[-zF(t)L 2 (t)(x 2 - x' 2 ) - iL(t)p • (x - x')]. (D.3) 

Note that because of the quadratic term in the exponent the integrations are nontrivial. 
For the two-particle density matrix we find analogously 

g^\x 1 ,x 2 ,x[,x' 2 ;t) =N(N-1) J dx 3 . . . dx N V*(x u x 2 , . . . , x N ; t) 

x ^(x[,x 2 , ...,x N ;t) 
1 (*) / xi x 2 x[ x' 2 \ 

= L(ty 92 \L(ty nty L(ty Lity ) 

x exp {-iF{t)(x\ + x\- xi - x 2 )) . (D.4) 
and the two-particle correlation function reads 

P?\x,y]t) = g { 2 \x,y,x,y;t) = jj^P^ (^'Z^)' ) ■ ( D ' 5 ) 

Other useful quantities such as non-equilibrium time-dependent correlation 
functions (e.g. n^'{p,t,t') ) or the multi-mode squeezing spectrum (S(k,k';t,t r ) = 
(n^(p^)n^ (p' \t')}) can also be easily computed using the scaling approach. 

Appendix E. Some technical details related to the derivation of ID and 2D 
momentum distribution at equilibrium 

Appendix E.l. Trapped weakly interacting Bose gases 

In order to describe a condensed Bose gas in a harmonic potential we adopt results 
of previous works [HI H3l H2] which consider phase fluctuations on top of the mean- 
field solution while density fluctuations are assumed to be negligible. This is a valid 
approximation for a sufficiently high number of weakly interacting particles at low 
temperatures. The temperature range where the density fluctuations are suppressed 
is T d ;§> T ;§> where the temperature of quantum degeneracy is T d = Nhu and 
T$ = T d hu /fx. 

Generically, the single-particle correlation can be represented as 

^(x,x') = y^x^xOexp (-i<(#x) - 0(x')) 2 >) , (E.l) 

where (0(x)) denotes the average over phase fluctuations. We assume the validity of 
the Thomas-Fermi approximation for the density 




where Rtf = y^p/^io/^o is the Thomas- Fermi radius. 
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In a ID geometry, taking into account thermal fluctuations and neglecting 
contributions from quantum fluctuations, one obtains the phase average [H] 



MX') - <Kx)n 



4T/i 



111 



Rtf ' 



Rtf - 



(E.3) 



Rt f Rt f ' 

For the 2D case an expression similar to the ID case can be derived. In this work 
we used the complete expression obtained by Xia et al. (Eq. (77) in Ref. [12]), which 
explicitly accounts for thermal and quantum fluctuations. As a result, at inter-particle 
distances much smaller than 2Rtf the correlations decay exponentially with a decay 
rate approximately given by rnksT/ (2irh 2 p(0)^ . However, for the dynamics studied in 
this paper we did not find significant effects from quantum corrections. 



Appendix E.2. One- and two-dimensional uniform Bose gases 

For a one-dimensional Bose gas it was recently shown [SB] that the effective field 
theory (Luttinger liquid) provides an extremely accurate description for a single-body 
correlation function at distances beyond the inter-particle separation. If we are not 
interested in its large momentum behavior it is legitimate to use this effective theory. 
The single particle correlation function in time-independent theory is then well known 
(see e.g. pQ). For nonzero temperatures it is given by (we omit oscillating terms) 

g?\x,x';0) = (&(x)$(x')) = Po B 

where £y — ftv s /T = h 2 irp/(m KT), p is the uniform equilibrium density, v s is the 
sound velocity, K is a Luttinger parameter which is related to the interaction strength 
c and B = {K/ti) 1 / 2K is Popov's factor. 

In the two-dimensional case, we consider a system below the Berezinskii-Kosterlitz- 
Thouless (BKT) transition. The correlation functions then decay algebraically with 
a temperature-dependent exponent, which tends to the universal value 1/4 when 
approaching the BKT transition from below. 



T 



po sinh(7r(x - x')/£ T ) 



(E.4) 



Appendix F. Dynamics of initially uniform systems 

Appendix F.l. Relating systems in the trap and without it 

The scaling approach can be used to establish a relationship between correlation 
functions in the model with time-dependent parameters (the system \1/) and the model 
with time-independent parameters (the system $). As we discussed in the main text, the 
trapping frequency Uo of the time-independent system is not fixed a priory. In particular, 
it can be put equal to zero from the very beginning. The scaling transformation therefore 
will relate the system in the time-dependent trap and a uniform system. The set of 
differential equations has to be modified accordingly. The aim of this appendix is to 
look into the behavior of the momentum distribution in this case. 
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The initial conditions state that the two wave functions are equal at t = 0. It means 
that the density distribution of the trapped system is homogeneous, corresponding to 
the uniform one. This is possible if we assume the existence of a length scale I on which 
this condition can be satisfied. Moreover we assume here that the Thomas-Fermi radius 
of a trapped system is large enough such that there is a finite region of x G [—l,l] D 
where the density is considered to be a constant. In the absence of a trapping potential 
this region is equal to the whole observation area. We assume that this region is large 
enough to contain a relatively large number of particles N. Using this length scale I as 
a sort of cut-off, we evaluate the momentum distribution in the finite window [— 1, 1} for 
examples of ID and 2D systems at finite temperature. 

Appendix F.2. Evaluation of momentum distributions in ID for the uniform system 

In the Luttinger liquid approximation at finite temperature we introduce £± = 
ir(x — d)/£t m terms of which g^(x,y;0) ~ (sinh£_) -1 / 2A ". This function decays 
exponentially at large distances and the limits of integration in £_ domain can be 
therefore extended from [— /, I] to (—oo, oo) to make analytic progress. One can 
easily realize that because of the additional structure in the exponent of Eq. (ID. II) . 



the expression for the momentum distribution is essentially different from the one at 
equilibrium. The corresponding integral is 

e -iO(m- I 



s*r(i-— ) (f.i) 



[sinhie-l] 1 / 2 ^ v 2K' 



x 



where C(t) = F(t)L 2 (t)^ + /n 2 + L(t)p£ T /n. The integration over is then performed 
in the finite interval [— /, /] corresponding to the size of the selected subsystem. The 
expression (IF. 21) is proportional to the equilibrium momentum distribution at t = 
provided that we take L(0) = 1. We find 

n(P,t) = 2 P0 L(t)^fr(l-±) 

w dc f n^-^(F(t)p(m + + P L(t))} 

-l/L(t) + ^ 

xT[l-^ - l -(F(t)L 2 (t)^+ P L(t))} 

_r[^ + i(F(t)L 2 m ++P L(t))} x 
r[i 

where L(t) = L(t)£ T /n. 

On the basis of this expression we have calculated a momentum distribution for 
various particular functions u(t) and m(t). Solving the set of consistency equations of 
Section [21 we obtained all the other functions v(t),L(t),F(t). This is illustrated in fig. 
IF II for particular choices of time-dependent functions u(t),m(t) and corresponds to a 



r[i- ± + UF(t)I?(t)t + + P L(t))] 



} (F.2) 
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4: 



~ 3 

Q_ 



— v s t/i; T =0.01 
•• v s t/^ T =0.5 
- v s t/^ T =0.75 
-v s t/^=1.0 




Figure Fl. The one-dimensional momentum distribution for different times computed 
within the Luttinger liquid model with K = 2, I = 10£l, normalized by a K- and po- 

dependent prefactor N Q = 2p (|f^) r(l — ). The trap frequency is changing 
linearly, ui(t) = t (j^-) , while the mass varies exponentially, m(t) = moe 2t " s ^ T . 



particular function v(t) found from solution of the Riccati equation. But additional 
simulations with various other choices of functions cu(t),m(t) suggest that the resulting 
momentum distribution defined as above in Eq. (I15j) has a step-like form. A formation of 
an effective momenta scale is associated with asymptotic emergence of microcanonical- 
type distribution. 

The Luttinger liquid expression for the gi-correlation function is a low-energy 
approximation for the true behavior of the correlation function. However, in the non- 
equilibrium dynamics we excite the whole spectrum and therefore the result for our time- 
dependent theory based on the exact equilibrium theory may appear to be different from 
the one based on the low-energy approximation. In what follows we demonstrate that 
the long-time behavior of the momentum distribution of the time-dependent system 
has a bounded support in momentum space. Our arguments can be applied to any 
exactly-solvable models. 

Suppose the g±(x, y)-correlation function is defined as a ground-state correlator 
of some field operators \E^(x): gi(x,y) = (x)^/(y)) . We also assume that 

the matrix elements of the operator ty(x) in the eigenbasis of the equilibrium problem 
are known. This implies that the form-factors F({A}, {/i}) = ({A}|\l/(0)|{ / u}) and the 
norms of the eigenstates |A) and are known. Here {/i}, {A} are the sets of numbers 
which characterize the eigenstates of a system of size 21. In particular these numbers 
can correspond to the solutions of the Bethe ansatz equations in the exactly-solvable 
problems. We also assume space- and time-translation invariance. Therefore the time- 
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IWI 2 IMI 2 



, (F.3) 



dependent gi function can be expanded as follows 

g{*\x, t- 0, 0) = £exp[z(P A - E„)t - z(P A - P M )a: 
M 

where Pa and P\ are, respectively, energy and momentum of the state |A). We assume 
also that the set {A} corresponds to the ground state. Introducing the coordinates 
£ = x — y and rj = x + y, the momentum distribution of the time-dependent system (we 
take for simplicity equal-time correlation function) can be written as 

fl/L(t) fl/L(t) 

i(p,t) = L(t)J2 / d£ / d V (FA) 



n 



04 



x exp 



l/L(t) 



i((P,-P x )+L(t)p + F(t)L 2 (t) V )£ 



l^(M,{A})p 



2 ||/i|| 2 



The ^-integration can be done easily, while after the ^-integration we obtain 
n(p,t) 



21 c ., 



iAii 2 ii/iii 2 



(F.5) 



where Si(z) is a sine-integral and x a = (P M — P\)l/L(t) +lp + aF(t)l 2 . The integrand is 
essentially proportional to F~ l (t) sin[(P M - P A +pL(t))/L(t)\/[{P ll - P\+ pL(t))/L(t)} 
and gives a main contribution to the sum when the momentum transfer is equal to 
pL(t). 

Appendix F. 3. Evaluation of the momentum distribution function in 2D for the 
uniform system. 

Here we evaluate the momentum distribution function for the 2D Bose gas below the 
BKT transition. We consider a system with time-dependent parameters and assume the 
validity of the long wavelength approximation. 



0.25 



0.15 



Q. 
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Figure F2. Normalized momentum distributions n(p) in the asymptotic-time limit 
for temperatures below the BKT-transition 
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According to (ID. 31) the momentum distribution in 2D is given by 

n(p,t) = L 2 (t) [ I drdr'g[*\r,r';0) (F.6) 

J\\r\\<l/L(t)J\\r'\\<l/L(t) 

x exp{-iF(t)L 2 (t)(r 2 - r' 2 ) - iL(t)p ■ (r - r')]. 

where the integration is restricted to a finite surface of the order of (2l/L(t)) 2 . We 
choose the density matrix in the scaling form corresponding to temperatures below the 
BKT transition, 

gi*\ r J-0) = PQ L^Y (F.7) 

where p is the density and rj = m T/ (2irh 2 p s (T)) (r] BKT = 1/4). Introducing the center 
of mass and relative coordinates 

Xi + X 2 V1+V2 {J?Q , 

x = x 1 -x 2 , y = yi-y2, x = — , Y= — - — , (F.8) 

and assuming the integration from —I to I we rewrite the momentum distribution as 

f l/L(t) f l/L(t) f l/L(t) f l/L(t) 

n(v,t)=L 2 (t)p v? dx dy dX dY (F.9 

J-l/L(t) J-UL(t) J-l/L(t) J-l/Ut) 



X 



-l/L(t) J-l/L(t) J-l/L(t) J-l/Ht) 

exp[i2F(t)L 2 (t)(xX + yY) + iL(t)(p x x + p y y)] 



(x 2 + y 2 )^ 2 

(p = (Px,Py))i which after integration over X and Y and changing variables to 
x -»■ x = 2lF(t)L(t)x = x, y ^y = 2lF{t)L{t)y = y, 



A = 2l 2 F(t), p x>y p x>y = (F.10) 
takes the following form 

4/ 2 Po^t f A , f A , sin(x) sm(y) e ip * x+ip y y 



nM = (2iF(t)ut))^ L dy[ 



(2lF(t)L(t)) 2 -v J-A J -a x y [x 2 + y 2 }n/ 2 

= E lap, (F.ll) 

a,/3=± 

l 2 p ^ f A , f A , sin[x(l + ap x )} sin[y(l + /%)] 
Q,/3 ~{2F{t)L{t)) 2 ^J-A J-A V x [x 2 + y 2 }v/ 2 y 

Now, using the integral / °° e _p2 "' M = p _1 / M r(l + -) we rewrite 



oo 



. ^^ N dt (F.12) 

and substitute back to Eq. (IF. 111) . Then the x and y integrals are separated now and 
can be performed using 

X) -^ = vrerf(J^Wn(C<) (F.13) 



x \2t x h t 

where we assume that the integration region can be effectively extended to infinity. This 
in particular is justified for large times when F(t) is a growing function of time or for 
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large I for arbitrary time. We therefore end up with the following integral 

n(P ' t} = (2lF(t)L(l)F%(l + *) a £ ± Sign(1 + ^ )Sign(1 + PPv) 

I ~, , . _ r /_ _ 1 

X 



/ erfdl + a^l^/^erfdl + ^l^ 1 /")^ (F.14) 

J 

fiables is transformed into the form 

/ x . r f°° eri(\a a \u)eri(\bB\u) , 

n(p,t) = J2 M«,f>J Q J+i ^ ( F - 15 ) 



which after the change of variables is transformed into the form 

erf(|a Q |M)erf(|6 j g|M 

where 

kt 7r 2 Z 2 po^(-7/)sign(a a )sign(fy) . , 

(2/F(t)L(t)) 2 -^r(l + f )2" 1 ' J 

a Q =l + ap x , b=l + f3f 3 p y (F.17) 

The last integral is equal to 



/(o*,^) = - ( |a Q r J B(--f,-,^) - V\b p \>>B(--£ 



2tt V 1 y a2'2'2 y 1 Pl v a 2 ' 2 '27 

+ M«2 (F , 8) 



where B(., .) is the Euler beta-function. So, finally we obtain 



n(p,t) 



(2lF(t)L(t)) 2 ^T(l + 2)2" 



2' 

x 5^ sign(l + ap x )sign(l + f3p y )I(a a , bp), (F.19) 

where J(a a , bp) is given in Eq. flF.18|> . a a = 1 + ap x /(2lF(t)) and 6 a = 1 + (3p y /(2lF(t)). 
We also introduced p XtV = p XtV /2lF(t). 

In fig. IF2I we plot the asymptotic behavior of the momentum distributions for 
various values of 77. Similarly to the one- dimensional case we find a step-like distribution 
which is smeared off when the BKT-transition is approached. 
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